1 Background

We compare the inference results from TMB, aghq, adam, and tmbstan. Import these inference results as follows:

tmb <- readRDS("depends/tmb.rds")
## Warning in checkMatrixPackageVersion(): Package version inconsistency detected.
## TMB was built with Matrix version 1.5.3
## Current Matrix version is 1.5.1
## Please re-install 'TMB' from source using install.packages('TMB', type = 'source') or ask CRAN for a binary version of 'TMB' matching CRAN's 'Matrix' package
aghq <- readRDS("depends/aghq.rds")
adam <- readRDS("depends/adam.rds")
tmbstan <- readRDS("depends/tmbstan.rds")

depends <- yaml::read_yaml("orderly.yml")$depends

Check that the parameters (latent field, hyperparameters, model outputs) sampled from each of the four methods are the same:

stopifnot(names(tmb$fit$sample) == names(aghq$quad$sample))
stopifnot(names(tmb$fit$sample) == names(adam$adam$sample))
stopifnot(names(tmb$fit$sample) == names(tmbstan$mcmc$sample))

1.1 Run details

For more information about the conditions under which these results were generated, see:

1.1.1 TMB

dependency_details <- function(i) {
  report_name <- names(depends[[i]])
  print(paste0("Inference results obtained from ", report_name, " with the query ", depends[[i]][[report_name]]$id))
  report_id <- orderly::orderly_search(query = depends[[i]][[report_name]]$id, report_name)
  print(paste0("Obtained report had ID ", report_id, " and was run with the following parameters:"))
  print(orderly::orderly_info(report_id, report_name)$parameters)
}

dependency_details(1)
## [1] "Inference results obtained from naomi-simple_fit with the query latest(parameter:tmb == TRUE)"
## [1] "Obtained report had ID 20230307-105805-9316189f and was run with the following parameters:"
## $tmb
## [1] TRUE
## 
## $aghq
## [1] FALSE
## 
## $k
## [1] 2
## 
## $ndConstruction
## [1] "product"
## 
## $tmbstan
## [1] FALSE
## 
## $niter
## [1] 1000
## 
## $nthin
## [1] 1
## 
## $adam
## [1] FALSE
## 
## $nsample
## [1] 5000
## 
## $area_level
## [1] 4

1.1.2 aghq

dependency_details(2)
## [1] "Inference results obtained from naomi-simple_fit with the query latest(parameter:aghq == TRUE && parameter:k == 3 && parameter:s == 8)"
## [1] "Obtained report had ID 20230413-123917-33a6956b and was run with the following parameters:"
## $aghq
## [1] TRUE
## 
## $k
## [1] 3
## 
## $s
## [1] 8
## 
## $grid_type
## [1] "pca"
## 
## $tmb
## [1] FALSE
## 
## $tmbstan
## [1] FALSE
## 
## $hmc_laplace
## [1] FALSE
## 
## $niter
## [1] 1000
## 
## $nthin
## [1] 1
## 
## $adam
## [1] FALSE
## 
## $nsample
## [1] 5000
## 
## $area_level
## [1] 4

1.1.3 adam

dependency_details(3)
## [1] "Inference results obtained from naomi-simple_fit with the query latest(parameter:adam == TRUE)"
## [1] "Obtained report had ID 20230322-225150-382d50aa and was run with the following parameters:"
## $adam
## [1] TRUE
## 
## $tmb
## [1] FALSE
## 
## $aghq
## [1] FALSE
## 
## $k
## [1] 2
## 
## $ndConstruction
## [1] "product"
## 
## $tmbstan
## [1] FALSE
## 
## $hmc_laplace
## [1] FALSE
## 
## $niter
## [1] 1000
## 
## $nthin
## [1] 1
## 
## $nsample
## [1] 5000
## 
## $area_level
## [1] 4

1.1.4 tmbstan

tmbstan_details <- dependency_details(4)
## [1] "Inference results obtained from naomi-simple_fit with the query latest(parameter:tmbstan == TRUE && parameter:niter > 50000)"
## [1] "Obtained report had ID 20230414-110613-89cd4244 and was run with the following parameters:"
## $tmbstan
## [1] TRUE
## 
## $niter
## [1] 1e+05
## 
## $nthin
## [1] 40
## 
## $tmb
## [1] FALSE
## 
## $aghq
## [1] FALSE
## 
## $k
## [1] 1
## 
## $grid_type
## [1] "product"
## 
## $s
## [1] 1
## 
## $hmc_laplace
## [1] FALSE
## 
## $adam
## [1] FALSE
## 
## $nsample
## [1] 5000
## 
## $area_level
## [1] 4
tmbstan_details
## $tmbstan
## [1] TRUE
## 
## $niter
## [1] 1e+05
## 
## $nthin
## [1] 40
## 
## $tmb
## [1] FALSE
## 
## $aghq
## [1] FALSE
## 
## $k
## [1] 1
## 
## $grid_type
## [1] "product"
## 
## $s
## [1] 1
## 
## $hmc_laplace
## [1] FALSE
## 
## $adam
## [1] FALSE
## 
## $nsample
## [1] 5000
## 
## $area_level
## [1] 4

1.2 Time taken

time_taken <- data.frame(
  "TMB" = tmb$time,
  "aghq" = aghq$time,
  "adam" = adam$time,
  "tmbstan" = tmbstan$time
)

write_csv(time_taken, "time_taken.csv")

time_taken
adam <- adam$adam

2 Histograms and ECDF difference plots

We create histograms and empirical cumulative distribution function (ECDF) difference plots of the samples from each method. All of the possible parameter names are as follows:

pars <- names(tmb$fit$sample)
pars
##  [1] "beta_alpha"                    "beta_anc_alpha"                "log_sigma_alpha_a"             "log_sigma_alpha_xa"            "u_alpha_xs"                    "log_sigma_rho_as"              "logit_phi_alpha_as"            "u_rho_xs"                      "u_rho_as"                     
## [10] "anc_clients_t1_out"            "population_t1_out"             "plhiv_attend_t1_out"           "anc_already_art_t1_out"        "anc_tested_pos_t1_out"         "anc_tested_neg_t1_out"         "ui_anc_rho_x"                  "log_sigma_ancrho_x"            "us_alpha_x"                   
## [19] "log_sigma_lambda_x"            "anc_rho_obs_t1_ll"             "hhs_prev_ll"                   "beta_anc_rho"                  "log_or_gamma"                  "u_alpha_a"                     "log_sigma_rho_xa"              "log_sigma_or_gamma"            "u_rho_xa"                     
## [28] "log_sigma_alpha_as"            "logit_phi_alpha_xs"            "logit_phi_rho_as"              "us_rho_xs"                     "anc_art_new_t1_out"            "untreated_plhiv_attend_t1_out" "untreated_plhiv_num_t1_out"    "artattend_ij_t1_out"           "log_betaT"                    
## [37] "anc_known_pos_t1_out"          "artattend_t1_out"              "OmegaT_raw"                    "ui_anc_alpha_x"                "ui_lambda_x"                   "log_sigma_ancalpha_x"          "u_rho_x"                       "u_alpha_x"                     "logit_phi_rho_x"              
## [46] "us_rho_x"                      "beta_rho"                      "u_alpha_xa"                    "logit_phi_alpha_a"             "u_rho_a"                       "logit_phi_rho_a"               "beta_lambda"                   "log_sigma_rho_a"               "u_alpha_as"                   
## [55] "us_alpha_xs"                   "log_sigma_alpha_xs"            "logit_phi_rho_xs"              "log_sigma_rho_xs"              "rho_t1_out"                    "anc_rho_t1_out"                "artnum_t1_out"                 "anc_alpha_t1_out"              "alpha_t1_out"                 
## [64] "infections_t1_out"             "anc_plhiv_t1_out"              "plhiv_t1_out"                  "lambda_t1_out"                 "log_sigma_alpha_x"             "logit_phi_alpha_x"             "log_sigma_rho_x"               "hhs_artcov_ll"                 "artnum_t1_ll"

We will produce plots about the following subset of them. There is no particular reason to choose this subset rather than other, it’s quite arbitrary.

pars_eval <- pars %in% c("beta_rho", "beta_alpha", "beta_lambda", "beta_anc_rho", "beta_anc_alpha", "u_rho_x", "us_rho_x", "u_rho_xs")
names(pars_eval) <- pars

We will write the tmbstan \(\hat R\) and ESS for the particular parameter on each plot to help in deciding how seriously to take the MCMC results:

rhats <- bayesplot::rhat(tmbstan$mcmc$stanfit)
ess_ratio <- bayesplot::neff_ratio(tmbstan$mcmc$stanfit)
niter <- 4 * tmbstan_details$niter / tmbstan_details$nthin
ess <- ess_ratio * niter

The beta_rho and beta_anc_rho parameters are used in other reports, so create dataframes here to save as artefacts:

beta_rho <- histogram_and_ecdf("beta_rho", i = 1, return_df = TRUE)
saveRDS(beta_rho$df, "beta_rho.rds")

beta_anc_rho <- histogram_and_ecdf("beta_anc_rho", return_df = TRUE)
saveRDS(beta_anc_rho$df, "beta_anc_rho.rds")

2.1 beta_rho

histogram_and_ecdf_helper <- function(par) lapply(1:sum(names(tmb$fit$obj$env$par) == par), histogram_and_ecdf, par = par)
histogram_and_ecdf_helper("beta_rho")
## [[1]]

## 
## [[2]]

2.2 beta_alpha

histogram_and_ecdf_helper("beta_alpha")
## [[1]]

## 
## [[2]]

2.3 beta_lambda

histogram_and_ecdf_helper("beta_lambda")
## [[1]]

## 
## [[2]]

2.4 beta_anc_rho

histogram_and_ecdf("beta_anc_rho")

2.5 beta_anc_alpha

histogram_and_ecdf("beta_anc_alpha")

2.6 logit

lapply(pars[stringr::str_starts(pars, "logit")], histogram_and_ecdf)

2.7 log_sigma

lapply(pars[stringr::str_starts(pars, "log_sigma")], histogram_and_ecdf)

2.8 u_rho_x

histogram_and_ecdf_helper("u_rho_x")
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

## 
## [[7]]

## 
## [[8]]

## 
## [[9]]

## 
## [[10]]

## 
## [[11]]

## 
## [[12]]

## 
## [[13]]

## 
## [[14]]

## 
## [[15]]

## 
## [[16]]

## 
## [[17]]

## 
## [[18]]

## 
## [[19]]

## 
## [[20]]

## 
## [[21]]

## 
## [[22]]

## 
## [[23]]

## 
## [[24]]

## 
## [[25]]

## 
## [[26]]

## 
## [[27]]

## 
## [[28]]

## 
## [[29]]

## 
## [[30]]

## 
## [[31]]

## 
## [[32]]

2.9 u_rho_x

histogram_and_ecdf_helper("u_rho_x")
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

## 
## [[7]]

## 
## [[8]]

## 
## [[9]]

## 
## [[10]]

## 
## [[11]]

## 
## [[12]]

## 
## [[13]]

## 
## [[14]]

## 
## [[15]]

## 
## [[16]]

## 
## [[17]]

## 
## [[18]]

## 
## [[19]]

## 
## [[20]]

## 
## [[21]]

## 
## [[22]]

## 
## [[23]]

## 
## [[24]]

## 
## [[25]]

## 
## [[26]]

## 
## [[27]]

## 
## [[28]]

## 
## [[29]]

## 
## [[30]]

## 
## [[31]]

## 
## [[32]]

2.10 us_rho_x

histogram_and_ecdf_helper("us_rho_x")
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

## 
## [[7]]

## 
## [[8]]

## 
## [[9]]

## 
## [[10]]

## 
## [[11]]

## 
## [[12]]

## 
## [[13]]

## 
## [[14]]

## 
## [[15]]

## 
## [[16]]

## 
## [[17]]

## 
## [[18]]

## 
## [[19]]

## 
## [[20]]

## 
## [[21]]

## 
## [[22]]

## 
## [[23]]

## 
## [[24]]

## 
## [[25]]

## 
## [[26]]

## 
## [[27]]

## 
## [[28]]

## 
## [[29]]

## 
## [[30]]

## 
## [[31]]

## 
## [[32]]

2.11 u_rho_xs

histogram_and_ecdf_helper("u_rho_xs")
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

## 
## [[7]]

## 
## [[8]]

## 
## [[9]]

## 
## [[10]]

## 
## [[11]]

## 
## [[12]]

## 
## [[13]]

## 
## [[14]]

## 
## [[15]]

## 
## [[16]]

## 
## [[17]]

## 
## [[18]]

## 
## [[19]]

## 
## [[20]]

## 
## [[21]]

## 
## [[22]]

## 
## [[23]]

## 
## [[24]]

## 
## [[25]]

## 
## [[26]]

## 
## [[27]]

## 
## [[28]]

## 
## [[29]]

## 
## [[30]]

## 
## [[31]]

## 
## [[32]]

2.12 us_rho_xs

histogram_and_ecdf_helper("us_rho_xs")

2.13 u_rho_as

histogram_and_ecdf_helper("u_rho_as")

2.14 u_alpha_x

histogram_and_ecdf_helper("u_alpha_x")

2.15 us_alpha_x

histogram_and_ecdf_helper("us_alpha_x")

2.16 u_alpha_xs

histogram_and_ecdf_helper("u_alpha_xs")

2.17 us_alpha_xs

histogram_and_ecdf_helper("us_alpha_xs")
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

## 
## [[7]]

## 
## [[8]]

## 
## [[9]]

## 
## [[10]]

## 
## [[11]]

## 
## [[12]]

## 
## [[13]]

## 
## [[14]]

## 
## [[15]]

## 
## [[16]]

## 
## [[17]]

## 
## [[18]]

## 
## [[19]]

## 
## [[20]]

## 
## [[21]]

## 
## [[22]]

## 
## [[23]]

## 
## [[24]]

## 
## [[25]]

## 
## [[26]]

## 
## [[27]]

## 
## [[28]]

## 
## [[29]]

## 
## [[30]]

## 
## [[31]]

## 
## [[32]]

2.18 u_alpha_a

histogram_and_ecdf_helper("u_alpha_a")

2.19 u_alpha_as

histogram_and_ecdf_helper("u_alpha_as")

2.20 u_alpha_xa

histogram_and_ecdf_helper("u_alpha_xa")

2.21 ui_lambda_x

histogram_and_ecdf_helper("ui_lambda_x")

2.22 ui_anc_rho_x

histogram_and_ecdf_helper("ui_anc_rho_x")

2.23 ui_anc_alpha_x

histogram_and_ecdf_helper("ui_anc_alpha_x")

2.24 log_or_gamma

histogram_and_ecdf_helper("log_or_gamma")

2.25 Outcome variables

These are all of the outcome variables which we might be interested in to monitor:

names(tmb$fit$sample)[!(names(tmb$fit$sample) %in% unique(names(tmb$fit$obj$env$par)))]
##  [1] "anc_clients_t1_out"            "population_t1_out"             "plhiv_attend_t1_out"           "anc_already_art_t1_out"        "anc_tested_pos_t1_out"         "anc_tested_neg_t1_out"         "anc_rho_obs_t1_ll"             "hhs_prev_ll"                   "u_rho_xa"                     
## [10] "anc_art_new_t1_out"            "untreated_plhiv_attend_t1_out" "untreated_plhiv_num_t1_out"    "artattend_ij_t1_out"           "anc_known_pos_t1_out"          "artattend_t1_out"              "rho_t1_out"                    "anc_rho_t1_out"                "artnum_t1_out"                
## [19] "anc_alpha_t1_out"              "alpha_t1_out"                  "infections_t1_out"             "anc_plhiv_t1_out"              "plhiv_t1_out"                  "lambda_t1_out"                 "hhs_artcov_ll"                 "artnum_t1_ll"

The most important ones (this has been checked with Jeff) are HIV prevalence (rho_t1_out), ART coverage (alpha_t1_out), and HIV incidence (lambda_t1_out).

  • HIV prevalence has 6417 variables (rows)
  • ART coverage has 6417 variables (rows)
  • HIV incidence has 6417 variables (rows)
dim(tmb$fit$sample$rho_t1_out)
## [1] 6417 5000
dim(tmb$fit$sample$alpha_t1_out)
## [1] 6417 5000
dim(tmb$fit$sample$lambda_t1_out)
## [1] 6417 5000

Within tmb$naomi_data there is a dataframe called mf_out (modelframe output) which contains the area, age, sex mapping these 6417 rows. As a sanity check, confirm that the number of areas times number of ages times number of sexes indeed equals the number of rows:

mf_out <- tmb$naomi_data$mf_out
(n_area <- length(unique(mf_out$area_id)))
## [1] 69
(n_age <- length(unique(mf_out$age_group)))
## [1] 31
(n_sex <- length(unique(mf_out$sex)))
## [1] 3
stopifnot(n_area * n_age * n_sex == dim(tmb$fit$sample$rho_t1_out)[1])

It makes sense only to assess inferential accuracy at the finest level so as to avoid double counting.

mf_out$area_id %>% unique()
##  [1] "MWI"           "MWI_1_1_demo"  "MWI_1_2_demo"  "MWI_1_3_demo"  "MWI_2_1_demo"  "MWI_2_2_demo"  "MWI_2_3_demo"  "MWI_2_4_demo"  "MWI_2_5_demo"  "MWI_3_1_demo"  "MWI_3_10_demo" "MWI_3_11_demo" "MWI_3_12_demo" "MWI_3_13_demo" "MWI_3_14_demo" "MWI_3_15_demo" "MWI_3_16_demo" "MWI_3_17_demo"
## [19] "MWI_3_18_demo" "MWI_3_19_demo" "MWI_3_2_demo"  "MWI_3_20_demo" "MWI_3_21_demo" "MWI_3_22_demo" "MWI_3_23_demo" "MWI_3_24_demo" "MWI_3_25_demo" "MWI_3_26_demo" "MWI_3_27_demo" "MWI_3_28_demo" "MWI_3_3_demo"  "MWI_3_4_demo"  "MWI_3_5_demo"  "MWI_3_6_demo"  "MWI_3_7_demo"  "MWI_3_8_demo" 
## [37] "MWI_3_9_demo"  "MWI_4_1_demo"  "MWI_4_10_demo" "MWI_4_11_demo" "MWI_4_12_demo" "MWI_4_13_demo" "MWI_4_14_demo" "MWI_4_15_demo" "MWI_4_16_demo" "MWI_4_17_demo" "MWI_4_18_demo" "MWI_4_19_demo" "MWI_4_2_demo"  "MWI_4_20_demo" "MWI_4_21_demo" "MWI_4_22_demo" "MWI_4_23_demo" "MWI_4_24_demo"
## [55] "MWI_4_25_demo" "MWI_4_26_demo" "MWI_4_27_demo" "MWI_4_28_demo" "MWI_4_29_demo" "MWI_4_3_demo"  "MWI_4_30_demo" "MWI_4_31_demo" "MWI_4_32_demo" "MWI_4_4_demo"  "MWI_4_5_demo"  "MWI_4_6_demo"  "MWI_4_7_demo"  "MWI_4_8_demo"  "MWI_4_9_demo"
mf_out$sex %>% unique()
## [1] "both"   "female" "male"
mf_out$age_group %>% unique()
##  [1] "Y000_004" "Y000_014" "Y000_064" "Y000_999" "Y005_009" "Y010_014" "Y010_019" "Y015_019" "Y015_024" "Y015_049" "Y015_064" "Y015_999" "Y020_024" "Y025_029" "Y025_034" "Y025_049" "Y030_034" "Y035_039" "Y035_049" "Y040_044" "Y045_049" "Y050_054" "Y050_064" "Y050_999" "Y055_059" "Y060_064"
## [27] "Y065_069" "Y065_999" "Y070_074" "Y075_079" "Y080_999"
mf_out_fine <- mf_out %>%
  tibble::rownames_to_column("id") %>%
  mutate(id = as.numeric(id)) %>%
  filter(
    area_id %in% paste0("MWI_4_", 1:32, "_demo"),
    sex %in% c("male", "female"),
    age_group %in%
      c(
        "Y000_004", "Y005_009", "Y010_014", "Y015_019", "Y020_024", "Y025_029",
        "Y025_034", "Y030_034", "Y035_039", "Y040_044", "Y045_049", "Y050_054",
        "Y055_059", "Y060_064", "Y065_069", "Y070_074", "Y075_079", "Y080_999"
      )
  )

An id column has been added to mf_out_fine to enable merging to the samples.

head(mf_out_fine)
nrow(mf_out_fine)
## [1] 1152
rho_t1_out_fine <- tmb$fit$sample$rho_t1_out[mf_out_fine$id, ]
dim(rho_t1_out_fine)
## [1] 1152 5000
rho_t1_out_fine[1:10, 1:10]
##              [,1]        [,2]        [,3]        [,4]        [,5]        [,6]        [,7]        [,8]        [,9]       [,10]
##  [1,] 0.003321178 0.003953769 0.003771124 0.003632820 0.003663655 0.003241938 0.002587445 0.003469123 0.002831204 0.003357822
##  [2,] 0.005144430 0.006124298 0.005841386 0.005627155 0.005674919 0.005021689 0.004007894 0.005373592 0.004385471 0.005201190
##  [3,] 0.006167904 0.007342715 0.007003518 0.006746667 0.006803933 0.006020744 0.004805257 0.006442658 0.005257952 0.006235956
##  [4,] 0.010151636 0.010907044 0.012744803 0.010384379 0.009977405 0.008354202 0.008791417 0.010521260 0.007694025 0.009703460
##  [5,] 0.025818257 0.027572174 0.026021053 0.025835757 0.024514982 0.024845195 0.018314273 0.025401685 0.021163558 0.023845284
##  [6,] 0.046337031 0.053570211 0.058109485 0.047331062 0.052884278 0.049100798 0.033744503 0.048900577 0.037903116 0.046274786
##  [7,] 0.057829702 0.068511303 0.068456822 0.062463838 0.066402044 0.059027290 0.044606395 0.062066116 0.048857612 0.062079671
##  [8,] 0.074018155 0.089557163 0.083031970 0.083779701 0.085443023 0.073009640 0.059906338 0.080610950 0.064287997 0.084342261
##  [9,] 0.097412216 0.126919810 0.107997937 0.107973605 0.117886716 0.107244600 0.081401087 0.105198444 0.089352703 0.103010505
## [10,] 0.112657742 0.133751227 0.118695841 0.127573009 0.118992506 0.100183649 0.090113330 0.108962172 0.090751596 0.103010068

Something weird happening with adam, go back and check this:

histogram_and_ecdf(par = "alpha_t1_out", i = 1)

3 KS plots

3.1 Individual parameters

ks_helper <- function(par, starts_with = FALSE, ...) {
  to_ks_df(par, starts_with = starts_with) %>% ks_plot(par, ...)
}

3.1.1 beta

ks_helper("beta", starts_with = TRUE)

3.1.2 logit

ks_helper("logit", starts_with = TRUE)
## Warning: Removed 116 rows containing missing values (`geom_tile()`).

3.1.3 log_sigma

ks_helper("log_sigma", starts_with = TRUE)
## Warning: Removed 116 rows containing missing values (`geom_tile()`).

3.1.4 u_rho_x

ks_helper("u_rho_x")

3.1.5 u_rho_xs

ks_helper("u_rho_xs")

3.1.6 us_rho_x

ks_helper("us_rho_x")

3.1.7 us_rho_xs

ks_helper("us_rho_xs")

3.1.8 u_rho_a

ks_helper("u_rho_a")

3.1.9 u_rho_as

ks_helper("u_rho_as")

3.1.10 u_alpha_x

ks_helper("u_alpha_x")

3.1.11 u_alpha_xs

ks_helper("u_alpha_xs")

3.1.12 us_alpha_x

ks_helper("us_alpha_x")

3.1.13 us_alpha_xs

ks_helper("us_alpha_xs")

3.1.14 u_alpha_a

ks_helper("u_alpha_a")

3.1.15 u_alpha_as

ks_helper("u_alpha_as")

3.1.16 u_alpha_xa

ks_helper("u_alpha_xa")

3.1.17 ui_anc_rho_x

ks_helper("ui_anc_rho_x")

3.1.18 ui_anc_alpha_x

ks_helper("ui_anc_alpha_x")

3.1.19 log_or_gamma

ks_helper("log_or_gamma")

3.1.20 Outcome variables

ks_df_out <- function(par) {
  to_ks_df(par = par, outputs = TRUE, id = mf_out_fine$id) %>%
  select(-par) %>%
  left_join(
    mf_out_fine %>%
      rename("full_id" = "id") %>%
      tibble::rowid_to_column("index") %>%
      mutate(index = as.numeric(index))
  )
}

ks_rho_t1_out <- ks_df_out(par = "rho_t1_out")
## Joining, by = "index"
ks_alpha_t1_out <- ks_df_out(par = "alpha_t1_out")
## Joining, by = "index"
ks_lambda_t1_out <- ks_df_out(par = "lambda_t1_out")
## Joining, by = "index"
ks_plot(ks_df = ks_rho_t1_out, par = "rho_t1_out", alpha = 0.2)

ks_plot(ks_df = ks_alpha_t1_out, par = "alpha_t1_out", alpha = 0.2)

ks_plot(ks_df = ks_lambda_t1_out, par = "lambda_t1_out", alpha = 0.2)
## Warning: Removed 116 rows containing missing values (`geom_tile()`).

3.2 Summary table

Values in dark blue are the lowest KS difference for that particular parameter:

options(dplyr.summarise.inform = FALSE)

ks_summary <- lapply(unique(names(tmb$fit$obj$env$par)), function(x) {
  to_ks_df(x) %>%
    group_by(method, index) %>%
    summarise(ks = mean(ks), par = x) %>%
    ungroup()
}) %>%
  bind_rows() %>%
  pivot_wider(names_from = "method", values_from = "ks") %>%
  rename(
    "Parameter" = "par",
    "KS(adam, tmbstan)" = "adam",
    "KS(aghq, tmbstan)" = "aghq",
    "KS(TMB, tmbstan)" = "TMB",
  )

r <- adam$quad$obj$env$random
x_names <- names(adam$quad$obj$env$par[r])
theta_names <- names(adam$quad$obj$env$par[-r])

dict <- data.frame(
  Parameter = c(unique(x_names), unique(theta_names)),
  Type = c(rep("Latent field", length(unique(x_names))), rep("Hyper", length(unique(theta_names))))
)

ks_summary <- ks_summary %>%
  left_join(dict, by = "Parameter")

ks_summary %>%
  gt::gt() %>%
  gt::fmt_number(
    columns = starts_with("KS"),
    decimals = 3
  ) %>%
  gt::tab_style(
    style = cell_fill(color = "#0072B2"),
    locations = cells_body(
      columns = `KS(adam, tmbstan)`,
      rows = `KS(adam, tmbstan)` < `KS(aghq, tmbstan)` & `KS(adam, tmbstan)` < `KS(TMB, tmbstan)`
    )
  ) %>%
    gt::tab_style(
    style = cell_fill(color = "#0072B2"),
    locations = cells_body(
      columns = `KS(aghq, tmbstan)`,
      rows = `KS(aghq, tmbstan)` < `KS(adam, tmbstan)` & `KS(aghq, tmbstan)` < `KS(TMB, tmbstan)`
    )
  ) %>%
    gt::tab_style(
    style = cell_fill(color = "#0072B2"),
    locations = cells_body(
      columns = `KS(TMB, tmbstan)`,
      rows = `KS(TMB, tmbstan)` < `KS(aghq, tmbstan)` & `KS(TMB, tmbstan)` < `KS(adam, tmbstan)`
    )
  )
index Parameter KS(adam, tmbstan) KS(aghq, tmbstan) KS(TMB, tmbstan) Type
1 beta_rho 0.077 0.045 0.061 Latent field
2 beta_rho 0.063 0.042 0.068 Latent field
1 beta_alpha 0.066 0.036 0.068 Latent field
2 beta_alpha 0.114 0.101 0.103 Latent field
1 beta_lambda 0.030 0.126 0.118 Latent field
2 beta_lambda 0.025 0.030 0.018 Latent field
1 beta_anc_rho 0.037 0.098 0.105 Latent field
1 beta_anc_alpha 0.042 0.070 0.059 Latent field
1 logit_phi_rho_x 0.562 0.562 0.134 Hyper
1 log_sigma_rho_x 0.642 0.642 0.157 Hyper
1 logit_phi_rho_xs 0.535 0.535 0.095 Hyper
1 log_sigma_rho_xs 0.669 0.669 0.277 Hyper
1 logit_phi_rho_a 0.508 0.508 0.043 Hyper
1 log_sigma_rho_a 0.566 0.566 0.103 Hyper
1 logit_phi_rho_as 0.545 0.545 0.092 Hyper
1 log_sigma_rho_as 0.570 0.570 0.126 Hyper
1 log_sigma_rho_xa 0.680 0.680 0.223 Hyper
1 u_rho_x 0.046 0.057 0.065 Latent field
2 u_rho_x 0.028 0.041 0.063 Latent field
3 u_rho_x 0.030 0.022 0.019 Latent field
4 u_rho_x 0.036 0.080 0.072 Latent field
5 u_rho_x 0.125 0.065 0.101 Latent field
6 u_rho_x 0.137 0.128 0.135 Latent field
7 u_rho_x 0.104 0.155 0.155 Latent field
8 u_rho_x 0.041 0.020 0.033 Latent field
9 u_rho_x 0.023 0.041 0.042 Latent field
10 u_rho_x 0.025 0.084 0.086 Latent field
11 u_rho_x 0.054 0.028 0.031 Latent field
12 u_rho_x 0.109 0.147 0.173 Latent field
13 u_rho_x 0.027 0.025 0.025 Latent field
14 u_rho_x 0.063 0.073 0.075 Latent field
15 u_rho_x 0.045 0.016 0.035 Latent field
16 u_rho_x 0.023 0.020 0.043 Latent field
17 u_rho_x 0.010 0.041 0.029 Latent field
18 u_rho_x 0.055 0.014 0.034 Latent field
19 u_rho_x 0.073 0.032 0.056 Latent field
20 u_rho_x 0.014 0.015 0.024 Latent field
21 u_rho_x 0.036 0.071 0.061 Latent field
22 u_rho_x 0.139 0.107 0.131 Latent field
23 u_rho_x 0.052 0.044 0.039 Latent field
24 u_rho_x 0.073 0.032 0.049 Latent field
25 u_rho_x 0.050 0.060 0.062 Latent field
26 u_rho_x 0.045 0.022 0.028 Latent field
27 u_rho_x 0.090 0.094 0.072 Latent field
28 u_rho_x 0.020 0.039 0.034 Latent field
29 u_rho_x 0.062 0.082 0.115 Latent field
30 u_rho_x 0.048 0.035 0.052 Latent field
31 u_rho_x 0.065 0.054 0.061 Latent field
32 u_rho_x 0.049 0.056 0.079 Latent field
1 us_rho_x 0.129 0.142 0.132 Latent field
2 us_rho_x 0.111 0.121 0.095 Latent field
3 us_rho_x 0.102 0.085 0.098 Latent field
4 us_rho_x 0.075 0.094 0.077 Latent field
5 us_rho_x 0.051 0.058 0.047 Latent field
6 us_rho_x 0.190 0.184 0.117 Latent field
7 us_rho_x 0.060 0.044 0.037 Latent field
8 us_rho_x 0.077 0.080 0.072 Latent field
9 us_rho_x 0.049 0.058 0.055 Latent field
10 us_rho_x 0.077 0.061 0.049 Latent field
11 us_rho_x 0.097 0.079 0.096 Latent field
12 us_rho_x 0.083 0.092 0.102 Latent field
13 us_rho_x 0.052 0.062 0.059 Latent field
14 us_rho_x 0.061 0.067 0.073 Latent field
15 us_rho_x 0.039 0.061 0.065 Latent field
16 us_rho_x 0.070 0.074 0.083 Latent field
17 us_rho_x 0.062 0.075 0.076 Latent field
18 us_rho_x 0.112 0.101 0.100 Latent field
19 us_rho_x 0.057 0.082 0.073 Latent field
20 us_rho_x 0.105 0.108 0.099 Latent field
21 us_rho_x 0.106 0.074 0.082 Latent field
22 us_rho_x 0.098 0.091 0.092 Latent field
23 us_rho_x 0.083 0.092 0.080 Latent field
24 us_rho_x 0.066 0.050 0.057 Latent field
25 us_rho_x 0.065 0.070 0.086 Latent field
26 us_rho_x 0.072 0.071 0.085 Latent field
27 us_rho_x 0.037 0.035 0.049 Latent field
28 us_rho_x 0.044 0.038 0.039 Latent field
29 us_rho_x 0.080 0.084 0.083 Latent field
30 us_rho_x 0.071 0.057 0.082 Latent field
31 us_rho_x 0.031 0.036 0.049 Latent field
32 us_rho_x 0.079 0.070 0.072 Latent field
1 u_rho_xs 0.088 0.081 0.090 Latent field
2 u_rho_xs 0.083 0.081 0.084 Latent field
3 u_rho_xs 0.113 0.098 0.126 Latent field
4 u_rho_xs 0.116 0.129 0.126 Latent field
5 u_rho_xs 0.215 0.160 0.180 Latent field
6 u_rho_xs 0.260 0.235 0.218 Latent field
7 u_rho_xs 0.189 0.178 0.186 Latent field
8 u_rho_xs 0.098 0.090 0.101 Latent field
9 u_rho_xs 0.119 0.112 0.141 Latent field
10 u_rho_xs 0.110 0.087 0.117 Latent field
11 u_rho_xs 0.108 0.089 0.094 Latent field
12 u_rho_xs 0.226 0.198 0.204 Latent field
13 u_rho_xs 0.109 0.109 0.110 Latent field
14 u_rho_xs 0.130 0.137 0.133 Latent field
15 u_rho_xs 0.157 0.131 0.159 Latent field
16 u_rho_xs 0.105 0.092 0.111 Latent field
17 u_rho_xs 0.096 0.083 0.088 Latent field
18 u_rho_xs 0.183 0.152 0.177 Latent field
19 u_rho_xs 0.161 0.147 0.167 Latent field
20 u_rho_xs 0.121 0.101 0.107 Latent field
21 u_rho_xs 0.116 0.132 0.121 Latent field
22 u_rho_xs 0.267 0.215 0.214 Latent field
23 u_rho_xs 0.169 0.143 0.156 Latent field
24 u_rho_xs 0.184 0.158 0.172 Latent field
25 u_rho_xs 0.114 0.092 0.113 Latent field
26 u_rho_xs 0.138 0.116 0.133 Latent field
27 u_rho_xs 0.186 0.182 0.187 Latent field
28 u_rho_xs 0.120 0.102 0.110 Latent field
29 u_rho_xs 0.209 0.171 0.201 Latent field
30 u_rho_xs 0.133 0.121 0.137 Latent field
31 u_rho_xs 0.152 0.133 0.154 Latent field
32 u_rho_xs 0.118 0.089 0.117 Latent field
1 us_rho_xs 0.030 0.031 0.031 Latent field
2 us_rho_xs 0.030 0.035 0.036 Latent field
3 us_rho_xs 0.022 0.022 0.026 Latent field
4 us_rho_xs 0.016 0.023 0.015 Latent field
5 us_rho_xs 0.026 0.022 0.032 Latent field
6 us_rho_xs 0.030 0.023 0.027 Latent field
7 us_rho_xs 0.038 0.042 0.048 Latent field
8 us_rho_xs 0.029 0.010 0.040 Latent field
9 us_rho_xs 0.033 0.014 0.030 Latent field
10 us_rho_xs 0.017 0.010 0.035 Latent field
11 us_rho_xs 0.036 0.016 0.042 Latent field
12 us_rho_xs 0.040 0.015 0.041 Latent field
13 us_rho_xs 0.019 0.023 0.031 Latent field
14 us_rho_xs 0.029 0.015 0.036 Latent field
15 us_rho_xs 0.031 0.021 0.039 Latent field
16 us_rho_xs 0.022 0.012 0.020 Latent field
17 us_rho_xs 0.020 0.022 0.019 Latent field
18 us_rho_xs 0.020 0.012 0.017 Latent field
19 us_rho_xs 0.025 0.028 0.010 Latent field
20 us_rho_xs 0.026 0.021 0.015 Latent field
21 us_rho_xs 0.016 0.024 0.014 Latent field
22 us_rho_xs 0.056 0.026 0.046 Latent field
23 us_rho_xs 0.041 0.025 0.031 Latent field
24 us_rho_xs 0.035 0.022 0.025 Latent field
25 us_rho_xs 0.053 0.042 0.053 Latent field
26 us_rho_xs 0.034 0.020 0.037 Latent field
27 us_rho_xs 0.024 0.019 0.024 Latent field
28 us_rho_xs 0.033 0.023 0.028 Latent field
29 us_rho_xs 0.033 0.023 0.028 Latent field
30 us_rho_xs 0.038 0.020 0.038 Latent field
31 us_rho_xs 0.036 0.017 0.038 Latent field
32 us_rho_xs 0.040 0.020 0.036 Latent field
1 u_rho_a 0.057 0.047 0.063 Latent field
2 u_rho_a 0.071 0.050 0.065 Latent field
3 u_rho_a 0.071 0.050 0.066 Latent field
4 u_rho_a 0.069 0.048 0.063 Latent field
5 u_rho_a 0.064 0.050 0.061 Latent field
6 u_rho_a 0.070 0.048 0.062 Latent field
7 u_rho_a 0.069 0.049 0.064 Latent field
8 u_rho_a 0.056 0.047 0.066 Latent field
9 u_rho_a 0.081 0.046 0.067 Latent field
10 u_rho_a 0.065 0.047 0.072 Latent field
1 u_rho_as 0.071 0.032 0.065 Latent field
2 u_rho_as 0.069 0.040 0.064 Latent field
3 u_rho_as 0.070 0.041 0.066 Latent field
4 u_rho_as 0.070 0.044 0.067 Latent field
5 u_rho_as 0.070 0.040 0.066 Latent field
6 u_rho_as 0.067 0.041 0.067 Latent field
7 u_rho_as 0.067 0.043 0.065 Latent field
8 u_rho_as 0.074 0.045 0.068 Latent field
9 u_rho_as 0.067 0.052 0.068 Latent field
10 u_rho_as 0.068 0.052 0.070 Latent field
1 logit_phi_alpha_x 0.509 0.509 0.100 Hyper
1 log_sigma_alpha_x 0.549 0.549 0.128 Hyper
1 logit_phi_alpha_xs 0.568 0.568 0.129 Hyper
1 log_sigma_alpha_xs 0.513 0.513 0.113 Hyper
1 logit_phi_alpha_a 0.553 0.553 0.101 Hyper
1 log_sigma_alpha_a 0.580 0.580 0.116 Hyper
1 logit_phi_alpha_as 0.537 0.537 0.072 Hyper
1 log_sigma_alpha_as 0.533 0.533 0.090 Hyper
1 log_sigma_alpha_xa 0.656 0.656 0.148 Hyper
1 u_alpha_x 0.125 0.116 0.147 Latent field
2 u_alpha_x 0.129 0.126 0.130 Latent field
3 u_alpha_x 0.136 0.112 0.119 Latent field
4 u_alpha_x 0.137 0.118 0.140 Latent field
5 u_alpha_x 0.130 0.098 0.128 Latent field
6 u_alpha_x 0.074 0.055 0.075 Latent field
7 u_alpha_x 0.073 0.060 0.061 Latent field
8 u_alpha_x 0.140 0.131 0.140 Latent field
9 u_alpha_x 0.164 0.144 0.159 Latent field
10 u_alpha_x 0.199 0.145 0.124 Latent field
11 u_alpha_x 0.147 0.126 0.149 Latent field
12 u_alpha_x 0.112 0.080 0.081 Latent field
13 u_alpha_x 0.138 0.104 0.124 Latent field
14 u_alpha_x 0.082 0.046 0.096 Latent field
15 u_alpha_x 0.115 0.079 0.089 Latent field
16 u_alpha_x 0.052 0.056 0.056 Latent field
17 u_alpha_x 0.085 0.116 0.121 Latent field
18 u_alpha_x 0.143 0.111 0.098 Latent field
19 u_alpha_x 0.109 0.101 0.126 Latent field
20 u_alpha_x 0.100 0.096 0.100 Latent field
21 u_alpha_x 0.170 0.145 0.103 Latent field
22 u_alpha_x 0.102 0.122 0.127 Latent field
23 u_alpha_x 0.121 0.103 0.116 Latent field
24 u_alpha_x 0.132 0.100 0.120 Latent field
25 u_alpha_x 0.082 0.108 0.115 Latent field
26 u_alpha_x 0.114 0.074 0.115 Latent field
27 u_alpha_x 0.123 0.114 0.130 Latent field
28 u_alpha_x 0.097 0.072 0.104 Latent field
29 u_alpha_x 0.086 0.082 0.126 Latent field
30 u_alpha_x 0.133 0.089 0.126 Latent field
31 u_alpha_x 0.081 0.063 0.089 Latent field
32 u_alpha_x 0.050 0.063 0.100 Latent field
1 us_alpha_x 0.064 0.066 0.067 Latent field
2 us_alpha_x 0.060 0.065 0.067 Latent field
3 us_alpha_x 0.080 0.071 0.068 Latent field
4 us_alpha_x 0.067 0.058 0.077 Latent field
5 us_alpha_x 0.103 0.062 0.105 Latent field
6 us_alpha_x 0.081 0.056 0.076 Latent field
7 us_alpha_x 0.023 0.014 0.032 Latent field
8 us_alpha_x 0.126 0.076 0.122 Latent field
9 us_alpha_x 0.116 0.091 0.118 Latent field
10 us_alpha_x 0.119 0.089 0.099 Latent field
11 us_alpha_x 0.099 0.086 0.086 Latent field
12 us_alpha_x 0.092 0.054 0.067 Latent field
13 us_alpha_x 0.103 0.075 0.100 Latent field
14 us_alpha_x 0.066 0.033 0.047 Latent field
15 us_alpha_x 0.090 0.061 0.080 Latent field
16 us_alpha_x 0.028 0.029 0.026 Latent field
17 us_alpha_x 0.056 0.064 0.081 Latent field
18 us_alpha_x 0.078 0.059 0.081 Latent field
19 us_alpha_x 0.074 0.070 0.097 Latent field
20 us_alpha_x 0.099 0.071 0.101 Latent field
21 us_alpha_x 0.118 0.083 0.101 Latent field
22 us_alpha_x 0.137 0.111 0.131 Latent field
23 us_alpha_x 0.079 0.063 0.071 Latent field
24 us_alpha_x 0.097 0.077 0.086 Latent field
25 us_alpha_x 0.081 0.073 0.089 Latent field
26 us_alpha_x 0.079 0.064 0.068 Latent field
27 us_alpha_x 0.070 0.069 0.068 Latent field
28 us_alpha_x 0.056 0.048 0.060 Latent field
29 us_alpha_x 0.105 0.082 0.096 Latent field
30 us_alpha_x 0.079 0.057 0.071 Latent field
31 us_alpha_x 0.053 0.037 0.044 Latent field
32 us_alpha_x 0.044 0.021 0.022 Latent field
1 u_alpha_xs 0.035 0.049 0.049 Latent field
2 u_alpha_xs 0.067 0.055 0.057 Latent field
3 u_alpha_xs 0.061 0.032 0.032 Latent field
4 u_alpha_xs 0.116 0.091 0.124 Latent field
5 u_alpha_xs 0.082 0.040 0.082 Latent field
6 u_alpha_xs 0.038 0.035 0.062 Latent field
7 u_alpha_xs 0.038 0.063 0.052 Latent field
8 u_alpha_xs 0.070 0.065 0.059 Latent field
9 u_alpha_xs 0.074 0.050 0.078 Latent field
10 u_alpha_xs 0.091 0.092 0.120 Latent field
11 u_alpha_xs 0.115 0.079 0.095 Latent field
12 u_alpha_xs 0.112 0.087 0.116 Latent field
13 u_alpha_xs 0.098 0.081 0.107 Latent field
14 u_alpha_xs 0.057 0.043 0.045 Latent field
15 u_alpha_xs 0.127 0.112 0.143 Latent field
16 u_alpha_xs 0.125 0.108 0.051 Latent field
17 u_alpha_xs 0.050 0.068 0.070 Latent field
18 u_alpha_xs 0.147 0.151 0.043 Latent field
19 u_alpha_xs 0.051 0.048 0.042 Latent field
20 u_alpha_xs 0.046 0.043 0.061 Latent field
21 u_alpha_xs 0.082 0.076 0.035 Latent field
22 u_alpha_xs 0.171 0.105 0.095 Latent field
23 u_alpha_xs 0.081 0.076 0.093 Latent field
24 u_alpha_xs 0.117 0.087 0.104 Latent field
25 u_alpha_xs 0.150 0.096 0.122 Latent field
26 u_alpha_xs 0.042 0.036 0.057 Latent field
27 u_alpha_xs 0.129 0.109 0.119 Latent field
28 u_alpha_xs 0.077 0.050 0.075 Latent field
29 u_alpha_xs 0.152 0.129 0.074 Latent field
30 u_alpha_xs 0.099 0.089 0.091 Latent field
31 u_alpha_xs 0.099 0.071 0.106 Latent field
32 u_alpha_xs 0.111 0.067 0.085 Latent field
1 us_alpha_xs 0.044 0.045 0.040 Latent field
2 us_alpha_xs 0.046 0.052 0.027 Latent field
3 us_alpha_xs 0.044 0.049 0.033 Latent field
4 us_alpha_xs 0.070 0.075 0.071 Latent field
5 us_alpha_xs 0.069 0.055 0.064 Latent field
6 us_alpha_xs 0.059 0.054 0.065 Latent field
7 us_alpha_xs 0.036 0.023 0.034 Latent field
8 us_alpha_xs 0.090 0.044 0.087 Latent field
9 us_alpha_xs 0.115 0.074 0.115 Latent field
10 us_alpha_xs 0.110 0.068 0.114 Latent field
11 us_alpha_xs 0.111 0.070 0.115 Latent field
12 us_alpha_xs 0.115 0.071 0.107 Latent field
13 us_alpha_xs 0.091 0.065 0.083 Latent field
14 us_alpha_xs 0.064 0.045 0.079 Latent field
15 us_alpha_xs 0.133 0.080 0.124 Latent field
16 us_alpha_xs 0.139 0.076 0.103 Latent field
17 us_alpha_xs 0.070 0.031 0.054 Latent field
18 us_alpha_xs 0.095 0.057 0.064 Latent field
19 us_alpha_xs 0.026 0.034 0.032 Latent field
20 us_alpha_xs 0.037 0.032 0.027 Latent field
21 us_alpha_xs 0.061 0.037 0.032 Latent field
22 us_alpha_xs 0.112 0.083 0.083 Latent field
23 us_alpha_xs 0.106 0.080 0.101 Latent field
24 us_alpha_xs 0.129 0.098 0.118 Latent field
25 us_alpha_xs 0.079 0.070 0.074 Latent field
26 us_alpha_xs 0.136 0.078 0.129 Latent field
27 us_alpha_xs 0.137 0.081 0.120 Latent field
28 us_alpha_xs 0.101 0.053 0.094 Latent field
29 us_alpha_xs 0.133 0.090 0.112 Latent field
30 us_alpha_xs 0.144 0.084 0.146 Latent field
31 us_alpha_xs 0.129 0.071 0.132 Latent field
32 us_alpha_xs 0.130 0.079 0.107 Latent field
1 u_alpha_a 0.059 0.038 0.062 Latent field
2 u_alpha_a 0.063 0.032 0.063 Latent field
3 u_alpha_a 0.061 0.029 0.064 Latent field
4 u_alpha_a 0.067 0.036 0.068 Latent field
5 u_alpha_a 0.064 0.035 0.068 Latent field
6 u_alpha_a 0.059 0.047 0.063 Latent field
7 u_alpha_a 0.064 0.043 0.062 Latent field
8 u_alpha_a 0.062 0.040 0.065 Latent field
9 u_alpha_a 0.063 0.033 0.063 Latent field
10 u_alpha_a 0.064 0.026 0.063 Latent field
11 u_alpha_a 0.063 0.035 0.059 Latent field
12 u_alpha_a 0.057 0.041 0.063 Latent field
13 u_alpha_a 0.067 0.035 0.061 Latent field
1 u_alpha_as 0.064 0.092 0.086 Latent field
2 u_alpha_as 0.099 0.097 0.104 Latent field
3 u_alpha_as 0.141 0.110 0.115 Latent field
4 u_alpha_as 0.133 0.109 0.113 Latent field
5 u_alpha_as 0.125 0.102 0.103 Latent field
6 u_alpha_as 0.107 0.087 0.091 Latent field
7 u_alpha_as 0.091 0.071 0.067 Latent field
8 u_alpha_as 0.059 0.061 0.056 Latent field
9 u_alpha_as 0.061 0.043 0.057 Latent field
10 u_alpha_as 0.062 0.044 0.063 Latent field
1 u_alpha_xa 0.049 0.032 0.039 Latent field
2 u_alpha_xa 0.032 0.023 0.025 Latent field
3 u_alpha_xa 0.016 0.017 0.017 Latent field
4 u_alpha_xa 0.055 0.051 0.067 Latent field
5 u_alpha_xa 0.106 0.068 0.109 Latent field
6 u_alpha_xa 0.100 0.089 0.090 Latent field
7 u_alpha_xa 0.057 0.063 0.058 Latent field
8 u_alpha_xa 0.060 0.055 0.055 Latent field
9 u_alpha_xa 0.042 0.050 0.046 Latent field
10 u_alpha_xa 0.111 0.114 0.117 Latent field
11 u_alpha_xa 0.112 0.098 0.104 Latent field
12 u_alpha_xa 0.108 0.092 0.100 Latent field
13 u_alpha_xa 0.060 0.058 0.064 Latent field
14 u_alpha_xa 0.071 0.083 0.044 Latent field
15 u_alpha_xa 0.033 0.042 0.050 Latent field
16 u_alpha_xa 0.044 0.053 0.035 Latent field
17 u_alpha_xa 0.032 0.021 0.029 Latent field
18 u_alpha_xa 0.025 0.036 0.010 Latent field
19 u_alpha_xa 0.095 0.063 0.082 Latent field
20 u_alpha_xa 0.024 0.030 0.040 Latent field
21 u_alpha_xa 0.089 0.090 0.067 Latent field
22 u_alpha_xa 0.227 0.166 0.183 Latent field
23 u_alpha_xa 0.032 0.021 0.019 Latent field
24 u_alpha_xa 0.044 0.021 0.049 Latent field
25 u_alpha_xa 0.099 0.055 0.067 Latent field
26 u_alpha_xa 0.024 0.018 0.019 Latent field
27 u_alpha_xa 0.101 0.086 0.087 Latent field
28 u_alpha_xa 0.045 0.026 0.051 Latent field
29 u_alpha_xa 0.045 0.037 0.021 Latent field
30 u_alpha_xa 0.051 0.051 0.049 Latent field
31 u_alpha_xa 0.015 0.027 0.013 Latent field
32 u_alpha_xa 0.025 0.021 0.034 Latent field
1 OmegaT_raw 0.512 0.512 0.018 Hyper
1 log_betaT 0.681 0.681 0.226 Hyper
1 log_sigma_lambda_x 0.630 0.630 0.164 Hyper
1 ui_lambda_x 0.110 0.106 0.096 Latent field
2 ui_lambda_x 0.106 0.103 0.107 Latent field
3 ui_lambda_x 0.107 0.104 0.102 Latent field
4 ui_lambda_x 0.102 0.100 0.104 Latent field
5 ui_lambda_x 0.106 0.102 0.122 Latent field
6 ui_lambda_x 0.089 0.095 0.112 Latent field
7 ui_lambda_x 0.124 0.108 0.108 Latent field
8 ui_lambda_x 0.102 0.102 0.093 Latent field
9 ui_lambda_x 0.110 0.104 0.109 Latent field
10 ui_lambda_x 0.159 0.164 0.130 Latent field
11 ui_lambda_x 0.121 0.129 0.121 Latent field
12 ui_lambda_x 0.092 0.099 0.113 Latent field
13 ui_lambda_x 0.125 0.129 0.124 Latent field
14 ui_lambda_x 0.106 0.105 0.100 Latent field
15 ui_lambda_x 0.116 0.114 0.112 Latent field
16 ui_lambda_x 0.120 0.104 0.122 Latent field
17 ui_lambda_x 0.126 0.108 0.109 Latent field
18 ui_lambda_x 0.094 0.105 0.106 Latent field
19 ui_lambda_x 0.117 0.098 0.119 Latent field
20 ui_lambda_x 0.126 0.114 0.099 Latent field
21 ui_lambda_x 0.113 0.103 0.096 Latent field
22 ui_lambda_x 0.098 0.097 0.108 Latent field
23 ui_lambda_x 0.100 0.094 0.095 Latent field
24 ui_lambda_x 0.103 0.103 0.100 Latent field
25 ui_lambda_x 0.105 0.108 0.107 Latent field
26 ui_lambda_x 0.131 0.133 0.131 Latent field
27 ui_lambda_x 0.099 0.104 0.101 Latent field
28 ui_lambda_x 0.093 0.101 0.094 Latent field
29 ui_lambda_x 0.097 0.106 0.095 Latent field
30 ui_lambda_x 0.101 0.097 0.099 Latent field
31 ui_lambda_x 0.103 0.120 0.096 Latent field
32 ui_lambda_x 0.109 0.109 0.112 Latent field
1 log_sigma_ancrho_x 0.502 0.502 0.016 Hyper
1 log_sigma_ancalpha_x 0.507 0.507 0.120 Hyper
1 ui_anc_rho_x 0.022 0.069 0.067 Latent field
2 ui_anc_rho_x 0.016 0.019 0.042 Latent field
3 ui_anc_rho_x 0.035 0.040 0.035 Latent field
4 ui_anc_rho_x 0.022 0.074 0.067 Latent field
5 ui_anc_rho_x 0.028 0.061 0.050 Latent field
6 ui_anc_rho_x 0.081 0.034 0.025 Latent field
7 ui_anc_rho_x 0.036 0.052 0.050 Latent field
8 ui_anc_rho_x 0.038 0.034 0.023 Latent field
9 ui_anc_rho_x 0.041 0.040 0.039 Latent field
10 ui_anc_rho_x 0.016 0.093 0.103 Latent field
11 ui_anc_rho_x 0.029 0.048 0.044 Latent field
12 ui_anc_rho_x 0.101 0.022 0.030 Latent field
13 ui_anc_rho_x 0.015 0.011 0.016 Latent field
14 ui_anc_rho_x 0.042 0.049 0.049 Latent field
15 ui_anc_rho_x 0.022 0.053 0.055 Latent field
16 ui_anc_rho_x 0.028 0.018 0.022 Latent field
17 ui_anc_rho_x 0.040 0.031 0.036 Latent field
18 ui_anc_rho_x 0.034 0.069 0.073 Latent field
19 ui_anc_rho_x 0.029 0.035 0.046 Latent field
20 ui_anc_rho_x 0.033 0.044 0.038 Latent field
21 ui_anc_rho_x 0.023 0.038 0.045 Latent field
22 ui_anc_rho_x 0.084 0.058 0.063 Latent field
23 ui_anc_rho_x 0.042 0.027 0.046 Latent field
24 ui_anc_rho_x 0.055 0.087 0.080 Latent field
25 ui_anc_rho_x 0.023 0.070 0.065 Latent field
26 ui_anc_rho_x 0.039 0.084 0.071 Latent field
27 ui_anc_rho_x 0.067 0.080 0.033 Latent field
28 ui_anc_rho_x 0.039 0.019 0.036 Latent field
29 ui_anc_rho_x 0.103 0.051 0.043 Latent field
30 ui_anc_rho_x 0.018 0.013 0.021 Latent field
31 ui_anc_rho_x 0.016 0.037 0.036 Latent field
32 ui_anc_rho_x 0.014 0.030 0.037 Latent field
1 ui_anc_alpha_x 0.042 0.047 0.068 Latent field
2 ui_anc_alpha_x 0.058 0.065 0.061 Latent field
3 ui_anc_alpha_x 0.099 0.082 0.085 Latent field
4 ui_anc_alpha_x 0.060 0.044 0.043 Latent field
5 ui_anc_alpha_x 0.042 0.063 0.053 Latent field
6 ui_anc_alpha_x 0.032 0.040 0.045 Latent field
7 ui_anc_alpha_x 0.062 0.076 0.056 Latent field
8 ui_anc_alpha_x 0.045 0.040 0.054 Latent field
9 ui_anc_alpha_x 0.068 0.068 0.075 Latent field
10 ui_anc_alpha_x 0.040 0.036 0.053 Latent field
11 ui_anc_alpha_x 0.039 0.041 0.034 Latent field
12 ui_anc_alpha_x 0.052 0.060 0.075 Latent field
13 ui_anc_alpha_x 0.035 0.029 0.038 Latent field
14 ui_anc_alpha_x 0.060 0.068 0.066 Latent field
15 ui_anc_alpha_x 0.076 0.074 0.092 Latent field
16 ui_anc_alpha_x 0.101 0.111 0.034 Latent field
17 ui_anc_alpha_x 0.099 0.119 0.121 Latent field
18 ui_anc_alpha_x 0.153 0.162 0.039 Latent field
19 ui_anc_alpha_x 0.086 0.082 0.067 Latent field
20 ui_anc_alpha_x 0.071 0.053 0.045 Latent field
21 ui_anc_alpha_x 0.076 0.058 0.056 Latent field
22 ui_anc_alpha_x 0.053 0.070 0.038 Latent field
23 ui_anc_alpha_x 0.047 0.056 0.065 Latent field
24 ui_anc_alpha_x 0.047 0.049 0.039 Latent field
25 ui_anc_alpha_x 0.101 0.066 0.081 Latent field
26 ui_anc_alpha_x 0.067 0.036 0.062 Latent field
27 ui_anc_alpha_x 0.063 0.079 0.064 Latent field
28 ui_anc_alpha_x 0.041 0.048 0.041 Latent field
29 ui_anc_alpha_x 0.161 0.175 0.071 Latent field
30 ui_anc_alpha_x 0.053 0.065 0.051 Latent field
31 ui_anc_alpha_x 0.050 0.046 0.064 Latent field
32 ui_anc_alpha_x 0.049 0.055 0.040 Latent field
1 log_sigma_or_gamma 0.518 0.518 0.020 Hyper
1 log_or_gamma 0.014 0.120 0.118 Latent field
2 log_or_gamma 0.012 0.060 0.071 Latent field
3 log_or_gamma 0.029 0.021 0.034 Latent field
4 log_or_gamma 0.025 0.048 0.050 Latent field
5 log_or_gamma 0.026 0.085 0.083 Latent field
6 log_or_gamma 0.015 0.065 0.052 Latent field
7 log_or_gamma 0.017 0.024 0.022 Latent field
8 log_or_gamma 0.040 0.022 0.027 Latent field
9 log_or_gamma 0.034 0.057 0.071 Latent field
10 log_or_gamma 0.024 0.022 0.029 Latent field
11 log_or_gamma 0.030 0.017 0.014 Latent field
12 log_or_gamma 0.023 0.102 0.105 Latent field
13 log_or_gamma 0.039 0.058 0.061 Latent field
14 log_or_gamma 0.026 0.019 0.011 Latent field
15 log_or_gamma 0.024 0.066 0.057 Latent field
16 log_or_gamma 0.027 0.049 0.043 Latent field
17 log_or_gamma 0.018 0.062 0.066 Latent field
18 log_or_gamma 0.034 0.077 0.087 Latent field
19 log_or_gamma 0.025 0.056 0.056 Latent field
20 log_or_gamma 0.036 0.025 0.035 Latent field
21 log_or_gamma 0.060 0.030 0.033 Latent field
22 log_or_gamma 0.021 0.024 0.021 Latent field
23 log_or_gamma 0.027 0.026 0.040 Latent field
24 log_or_gamma 0.024 0.051 0.057 Latent field
25 log_or_gamma 0.027 0.095 0.081 Latent field
26 log_or_gamma 0.013 0.013 0.021 Latent field
27 log_or_gamma 0.021 0.013 0.018 Latent field
28 log_or_gamma 0.020 0.027 0.025 Latent field
29 log_or_gamma 0.024 0.041 0.031 Latent field
30 log_or_gamma 0.028 0.056 0.040 Latent field
31 log_or_gamma 0.014 0.053 0.052 Latent field
32 log_or_gamma 0.021 0.052 0.039 Latent field
 
ks_summary %>%
  filter(Type != "Hyper") %>%
  ks_plot_many("TMB", "aghq")
## Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as of ggplot2 3.3.4.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.
## Warning: Removed 116 rows containing missing values (`geom_tile()`).
## Picking joint bandwidth of 0.00411

(ks_summary_out <- ks_summary %>%
  group_by(Type) %>%
  summarise(
    TMB = mean(`KS(TMB, tmbstan)`),
    aghq = mean(`KS(aghq, tmbstan)`),
    adam = mean(`KS(adam, tmbstan)`)
  )) %>%
  gt::gt()
Type TMB aghq adam
Hyper 0.11711667 0.56767500 0.56767500
Latent field 0.07540086 0.06824325 0.07582955

Save some things for output:

saveRDS(ks_summary_out, "ks_summary_out.rds")

3.3 Investigation into large KS values

Want to create rank order lists of largest KS differences between methods:

ks_comparison <- function(ks_summary, method1, method2, n) {
  ks_method1 <- paste0("KS(", method1, ", tmbstan)")
  ks_method2 <- paste0("KS(", method2, ", tmbstan)")
  ks_summary[["KS difference"]] <- ks_summary[[ks_method1]] - ks_summary[[ks_method2]]
  ks_summary %>%
    filter(Type != "Hyper") %>%
    arrange(desc(`KS difference`)) %>%
    head(n = n)
}

Nodes where TMB beats aghq:

ks_comparison(ks_summary, method1 = "aghq", method2 = "TMB", 10)
histogram_and_ecdf(par = "ui_anc_alpha_x", i = 18)

histogram_and_ecdf(par = "u_alpha_xs", i = 18)

histogram_and_ecdf(par = "ui_anc_alpha_x", i = 29)

Nodes where aghq beats TMB:

ks_comparison(ks_summary, method1 = "TMB", method2 = "aghq", 10)  
histogram_and_ecdf(par = "us_alpha_xs", i = 30)

histogram_and_ecdf(par = "us_alpha_xs", i = 31)

histogram_and_ecdf(par = "us_alpha_xs", i = 26)

3.4 Correlation between KS values and ESS

Is there any correlation between the value of \(\text{KS}(\texttt{method}, \texttt{tmbstan})\) for a particular parameter and the ESS of that parameter from tmbstan output?

ks_summary %>%
  pivot_longer(
    col = starts_with("KS"),
    names_to = "Method",
    values_to = "KS"
  ) %>%
  mutate(
    Method = fct_recode(Method,
      "adam" = "KS(adam, tmbstan)",
      "aghq" = "KS(aghq, tmbstan)",
      "TMB" = "KS(TMB, tmbstan)"
    )
  ) %>%
  group_by(Parameter) %>%
  mutate(
    par_num = case_when(
      max(index) > 1 ~ paste0(Parameter, "[", index, "]"),
      TRUE ~ Parameter
    )
  ) %>%
  ungroup() %>%
  left_join(data.frame(ess) %>%
    tibble::rownames_to_column("par_num")
  ) %>%
  ggplot(aes(x = ess, y = KS)) +
    geom_point(alpha = 0.5) +
    geom_smooth(method = "lm", color = "#56B4E9") +
    facet_grid(~ Method) +
    theme_minimal() +
    labs(x = "ESS", y = "KS(method, tmbstan)")
## Joining, by = "par_num"
## `geom_smooth()` using formula = 'y ~ x'

4 MMD

#' To write!

5 PSIS

Suppose we have two sets of samples from the posterior. For each sample we are going to want to evaluate the log-likelihood, so that we can calculate the log-likelihood ratios. We can extract the TMB objective function for the log-likelihood as follows:

tmb$fit$obj$fn()
## [1] NaN
#' To write!